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Abstract. We study the morphological evolution of surfaces during ion sputtering 
and we compare their dynamical roughening with aeolian ripple formation in sandy 
deserts. We show that, although the two phenomena are physically different, they 
must obey to similar geometrical constraints and therefore can be described within 
the same theoretical framework. The present theory distinguish between atoms that 
stay bounded in the bulk and others that are mobile on the surface. We describe 
the excavation mechanisms, the adsorption and the surface mobility by means of a 
continuous equation derived from the study of dune formation on sand. This approach 
can explain the different dynamical behaviors experimentally observed in metals or in 
semiconductors and amorphous systems. We also show that this novel approach can 
describe the occurrence of ripple rotation in the {x, y) plane induced by changes in the 
sputtering incidence angle. 



PACS numbers: 45.70.-n 61.43.-j 68.35.Ct 68.55.Jk 




Figure 1. Ripples on sand (Gobi desert) and Ripples on surfaces (Ag under ion 
sputtering). The line is indicative of the different length scales involved: it represents 
1 m in the left picture and 50 nm in the right one. 



1. Introduction 

When an ion liits a surface liberates locally a large amount of energy that melts a region 
of the solid immediately below. For geometrical reasons, the sputtering effect depends 
on the surface- curvature: the energy concentrates on regions of positive curvature and 
this favorites the excavation of valleys and the growth of hills. On the other hand, 
thermal diffusion and surface tension tend to smoothen the irregularities by flattening 
the surface. It has been observed that under the combined action of these mechanisms 
the surface tends to create spontaneously ripples [HllllSllllElinilZllHl- In nature, ripples 
are commonly observed in sandy-deserts as the result of dynamical instability of the sand 
surface under the action of a sufficiently strong wind In this case, the formation of 
ripples is commonly associated with the effect produced by some grains that are lifted 
from the sand-bed and accelerated by the wind. These grains, when re-impact with the 
bed, splash up a number of other grains. Most of these grains return to the bed leading to 
a local rearrangement, whereas some other are accelerated by the wind and impact again 
on the bed after a certain 'saltation' length. In the literature, many studies have been 
devoted to understanding the mechanism of ripple formation fOl El El El El ■ In 
particular, an hydrodynamical model for aeolian ripple formation, based on a continuum 
dynamical description with two species of grains (immobile and rolling grains), was 
proposed with success by Bouchaud et al. jTHl El El E] • The main ingredient of such 
a model is a bilinear differential equation, for the population of the two species of grains, 
which shows the instability of a flat bed against ripple formation. 

In this paper we show that the same reasoning which has been used to describe the 
sand ripples formation in deserts, applied to the studies of dynamical surface roughening. 
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leads to an accurate description of the morphogenesis and evolution of ripples on crystal 
and amorphous surfaces during ion sputtering. The present approach contains the 
Bradley-Harper approach [2111 and its nonlinear extension by Rost et al. j2I] and by Park 
et al. j22j (based on a Kuramoto-Sivashisky and Kardar-Parisi-Zhang type equations 
jSniElI) and it is also able to describe some of the crucial experimental features observed 
in these systems 0. In particular, by means of this approach we can describe the two 
distinct dynamical behaviors experimentally observed in amorphous/semiconductors 
systems and in metals In the first case (amorphous/semiconductors) we find that the 
ripples growth exponentially fast at constant wavelength A up to a critical roughening 
Wc at which the growing process interrupts. On the other hand, in metals (when 
the Erlich-Schwoebel barrier is active), we can observe a transition between an initial 
exponential to a slower power-law growth of the roughness (the root mean square of the 
height profile) . In this regime the ripple-wavelength tends also to growth with time and 
re-nucleation effects are observed. Finally, we show that the intriguing phenomena of 
"ripple rotation" [0] , associated with a change in the sputtering incidence angle, can be 
conveniently described by means of the present approach. 

The main novelty of this work is the use of a new kind of equation to describe 
the surface instabilities under ion sputtering. In the present approach we clearly 
distinguish between three mechanisms: 1) Erosion; 2) Adsorption-Condensation; 3) 
Mobility. Physically, all these three mechanisms take place simultaneously at a surface 
under ion-bombardment. The possibility to distinguishing between them clarifies the 
approach, simplifies the choice of the parameters and enlarges the descriptive power of 
the theory. Moreover, we distinguish between atoms that stay in the bulk in bounded 
position and others that stay on the surface and are mobile. This (simplified) view is 
more realistic than considering the whole system as a fluid. We note that these novel 
elements of the present theory do not introduce -per se- extra parameters in the analogy 
with the well established Bradley-Harper theory [20! • However, our choice in this paper 
has been to study the equations in their most complete form introducing therefore extra 
parameters. We show that, from the linear analysis, the Bradley-Harper solution for 
the ripple wavelength as well as other solutions for typical ripple periodicity in sandy 
deserts can be all retrieved and extended. 

2. Particle mobility and Ripple dynamics 

When the surface of a solid is taken under ion sputtering some atoms in the proximity of 
the surface receive energy from the sputtered ions and pass from a bounded - 'immobile' 
- solid state to a 'mobile' melted state. The opposite mechanism is also allowed: some 
mobile atoms can gain in energy by becoming immobile and bounding in a given position 
in the solid. A certain fraction of atoms might also be dispersed into the atmosphere. 
Let us call h{r, t) the height of surface profile made of immobile -bounded- atoms and 
call -R(r, t) the height of mobile -melted- atoms. In analogy with the theory developed 
to explain the dynamical evolutions of dunes in deserts [TBI El UHl CHI , we describe the 
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mechanisms of excavation, exchange between mobile and immobile atoms and surface 
displacement of mobile atoms in term of the following differential equation: 
dh 

— = -T{R,h)ex + T{R,h)ad 

f)R 

— = - VJ(i?, h) + {l- <j))T{R, h),, - T{R, h)ad . 

(1) 

Where T{R,h)ex is the rate of atoms that are excavated under the action of the 
sputtering, and (1 — cf)) is the part of them that pass from immobile to mobile, whereas 
(j) is the fraction that is dispersed into the atmosphere. 

Let us now write in details the various terms contained in EqC] 



2.1. Excavation 

The excavation effect must clearly depend on the number and velocity of the sputtered 
ions (i.e. its flux), but also the local shape and orientation of the surface might play 
an important role. Indeed, the energy transmitted by the impacting ions concentrate 
more in regions of the surface with positive curvature. Moreover, part of the surface 
facing the flux are likely to experience a different erosion respect to others which are less 
exposed to the flux. Such dependence of the sputtering yields from the combination of 
the surface curvature and the incoming ion flux orientation has been described in great 
details in the literature (see for instance ^ and references therein). Here we adopt 
the same reasoning that lead to the Bradley and Harper equation [20] • Crystalline 
orientation and anisotropies might be also taken into account. We can write: 

r(i?,/.)e. = r^(l + aV/. + (6.^ + 6,^)/.) . (2) 

Here r/ is the sputtering flux; a vectorial parameter associated with the 

flux-direction-dependent erosion and 6^, hy are associated with the curvature-dependent 
sputtering erosion. All these term are implicitly dependent on the orientation of the 
incoming ion flux, this in analogy with PU] . 

A noise term ^{x,y,t) can be eventually added to Eq|2] This term might mimic 
the fact that the incoming flux is made -after all- of discrete particles and that in 
experiments certain amounts of noise and randomness are unavoidable. 



2.2. Adsorption 

The rate of adsorption of mobile atoms into immobile solid positions must be dependent 
on the quantity of mobile atoms in a given spatial position. Similarly to the excavation 
process, the adsorption is also dependent on the local curvature and orientation. We 
can write: 

T{R,h)ad = R{l + c\/h+{dx— + dy—)hJ , (3) 
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where the parameter 7 is the recombination rate and c = {cx,Cy) and dx, dy are 
associated to the different probabihties of recombination in relation with the local 
orientation and shape of the surface. 

Note that Eqs. El El contain the same terms [23] as the ones proposed in the 
literature for the formation of aeolian dunes in the so-called hydrodynamical model 
[ini nil im nil I2ni- indeed, in deserts, sand grains are lifted from the sand-bed 
and readsorbed into it with a probability which is dependent on the local shape and 
orientation of the dunes. Eqs. El represent the simplest analytical expressions which 
formally take into account these shape and orientation dependences. In the quest for 
simple explanations, such equations are therefore rather universal. 



2.3. Mobility 

Mobile atoms will move on the surface, and the quantity J(-R, h) in EqQis the 'current' 
of these atoms. In surface growth, there are two main mechanisms that are commonly 
indicated as responsible for the surface mobility of atoms [2ZI. The first is a current, 
driven by the variations of the local chemical potential, which tends to smoothen the 
surface asperity moving atoms from hills to valleys. The second is a current induced by 
the Erlich-Schwoebel barrier which -on the contrary- moves atoms uphill. In addiction 
to these main mechanisms we might also have to take into account a drift velocity and 
a random thermal diffusion, obtaining: 

V/) 

J{R,h) = ICRV(y^h) + sR -—— + ^rR-VVR . (4) 

1 + [adVhy 

In this equation, the first term describes a deterministic diffusion driven by the variations 
of the chemical potential which depends on the local shape of the surface; the second 
term is associated with the uphill current due to the Erlich-Schwoebel barrier and 
is a constant associated with the characteristic length. The quantity v = {vx.Vy) is a 
drift velocity of the mobile atoms on the surface, whereas T) is the dispersion constant 
associated with the random thermal motion. (The coefficient T) is related with a non- 
deterministic diffusion mechanism and it could mimic the evaporation-condensation 
effects.) 

Note that EqEl is substantially different from the one proposed in the literature to 
describe ripples in granular media [121 El UHl I2H1 121] ■ Here the current is supposed to 
be dependent on the local shape and orientation of the surface (the h{r,t) profile). The 
equations describing sandy deserts can be retrieved from Eq|3]by imposing /C = and 
s = 0, but -on the contrary- in surface growth these two parameters are the leading 
terms of the equation and play the role of control parameters in the dynamics of ripple 
formation. Nonetheless, these terms describe a rather simple dependence of the dynamic 
of particles on a surface on the geometrical shape of the surface itself. Again, in our 
seek for universality, we expect that similar terms can be profitably introduced in the 
context of aeolian sand ripples in order to describe specific phenomena (associated, for 
instance, with packing properties [30] or granular flow [M]) which relate the current of 
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Figure 2. Numerical solutions of Eq^at various times |S2I indicate that under the 
action of ion sputtering the surface develops an instability which leads to the formation 
of ripples with a well defined characteristic wavelength. In the figure the black-tick 
line is the final surface-profile, whereas the thinner gray lines (green online) are some 
profiles at previous times. See Appendix [Appendix D| for details. 
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Figure 3. A numerical solution of Eq^in two dimensions. See Appendix |Appendix D| 
for details. 

grains witli tlie dune-shapes. 

It should be noted that the factors a, c and v in Eqs l2f HI and El are vectors (i.e. they 
have two distinct components in the {x,y) directions). Indeed, off- normal sputtering 
introduce anisotropy, and crystal surfaces are in general anisotropic. Therefore one 
must take into account the dependence of the parameters on the relative orientation 
of both the crystal-surface and the sputtering direction. In EqE] the deterministic 
diffusion parameter /C is assumed isotropic. On the other hand, on crystalline 
surfaces anisotropics are expected and sometimes might play an important role |22j . 
The extension of the results presented in this paper to the case of anisotropic /C is 
straightforward. 
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3. Dispersion Relation 

A trivial solution of Eq^ can be written for a completely flat surface: h{r,t) = ho{t) 
and R{r,t) = Rq. In this case, we obtain Rq = {1 — (p)ri/'y and ho{t) = —(prjt + const.. 
This describes a surface that rests flat and it is eroded with a speed equal to ^r^. 
But this behavior is only hypothetical since -in general- the dynamics of the surface- 
profile presents instabilities against spontaneous roughening and therefore its evolution 
is more complex. For instance, numerical solutions of Eq^ are shown in FigEl (for the 
1-dimensional case) and in FigsOl El IA2I (for the two-dimensional case). We observe 
that, in a certain range of the parameters, the surface is unstable and periodic ripples 
or other instabilities are formed spontaneously. 

3.1. Stability analysis 

In order to infer indications about the amplification or the smoothing of small 
perturbations and to deduce an analytical expression for the ripples wave-length at 
their beginning, we performe a stability analysis on Eq^ For this purpose we assume 
that the surface-profile is made by the combination of a fiat term plus a rough part: 

R{r,t) = Ro + Ri{r,t) 

h{r,t) = ho{t) + h{r,t) , (5) 

with -Ri(r, t) = Ri exp^itut + ikr) and hi(v, t) = hi exp{iujt + ihr). We substitute these 
quantities into Eq^ and linearize the equation by taking only the first order in Ri and 
hi. A Fourier analysis (see Appendix [Appendix A| ) shows that such a linearized equation 
admits solutions when the frequencies u and the wave vectors k satisfy: 

[icj + 7 + ikv + k"^!)] ■ 

{tUJ + tk [Vi - (1 - 0)V2] - [Vi - (1 - <P)V2]} - 

7(1 - 0) [«k(vi-V2)-A;2 {V1-V2 - Si)-k^JCi] = ; 

(6) 

where, to simplify the equations, we have introduced the following notation: 

vi = r/a si = 775/7 fCi = 77/0/7 

T^i = v{bx cos(a)^ + by sin(Q;)^) 
V2 = ^{d^cos{af + dysm{af) . 

Where, a is the azimuthal angle corresponding the given direction of k in the (x, y) plane 
(therefore kx = |k| cos(q;), ky = |k| sin(a)). Equation IHl establishes a dispersion relation 
a;(k) that is a complex function with two branches corresponding to the solutions of the 
quadratic Equation El 

The dispersion relation EqjH] reduces to the one from the linear analysis in the 
Bradley-Harper theory jSOjj when the coefficients V, v, s, ay, 7, c, d^, dy, si, V2 and V2 
are set equal to zero. 
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Figure 4. The imaginary part of tlie dispersion relation Im{uj) can assume 
negative values which are associated with the surface instability (arbitrary units) . The 
amplitude of modes with wavelengths A > 27r/fc* will grow exponentially fast. The tick 
line is the imaginary part of the analytical solution of Eq[S| whereas the thinny-gray 
line is the approximated expression (at the fourth order in k) obtained for small ion 
flux (7/ small). (To draw this figure we set: 7 = 8, </> = 0.5, v = 0.4, V = 0.2, ICi = 0.06, 
si = 0.09, vi = 0.01, V2 = 0.001, Vi = 0.015 and 2?2 = 0.023.) 



4. Surface Instabilities 



The kinetic growth of the surface instabihty is related to the immaginary part of co'(k). 
Indeed, /m(a;(k)) corresponds to modes with amplitudes that change exponentially fast 
in time, and negative values correspond to unstable modes that increase with the time. 
We can therefore study Im{uj) from the solution of EqlHland search for the region of k 
in which Im{uj) is negative. The most unstable mode is the one that grows faster and 
it corresponds to the value of k at which Im{u}) reaches its most negative value (see 
FigsHandEH). 

The solution of EqlHlfor Im{(jj), is 



2Im{(jj)± = 7 + 



± 



Ai + (Af + 4Ai)V2 



(7) 



where we have 



Ai = y - 



Vl 



jV2 



27 



P - (1 - 20)Pi + (1 - 0)^2 + 2(1 - k 



I [(P + Pi - (1 - 0)1)2] ' - 47(1 - 0)/Ci}fc 
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Figure 5. Surface-plot and contour lines plot for the imaginary part of the dispersion 
relation Im{ijj) from the analytical solution of Eq^lin the two-dimesional case. Figures 
(a) and (b) show the occurrence of a 90 degrees rotation in the orientation of the 
instability mode associated with the relative variations of the coefficients h,j, and hy. 
To obtain the plots on the left side we set: v, vi, V2 and s = 0, 7 = 0.05, 77 = 1, 



V = 0.1, b,^ = 0.25, by = 0.5, d^; = 0.025, dy = 0.025, /Ci = 0.833, 
plots uses the same set of parameters with bx and by exchanged. 



0.2. The right 



and 



(8) 



A, 



7 



V + (1 - 20)Vi - (1 - 0)V2 



+ 



I? + Pi - (1 - <P)V2 



V - Vi + (1 



JV2 



(9) 

Let us first observe that in absence of sputtering (i.e. when = and therefore, 
vi = 0, V2 = 0, Vi = 0, V2 = 0, Si = 0, /Ci = 0) the solutions of EqElare u;(k) = 
and = — kv + 2(7 + k'^V). In this case, the imaginary part of uj{k) is non- 

negative, therefore we -correctly- expect no spontaneous corrugation of the surface. 
On the contrary, when the sputtering is active (77 7^ 0), the immaginary part of uj{k.) 
can assume negative values. This is shown in FigE] where a plot of Im{uj)^ is reported 
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(along a given direction of the vector k). As one can see in FigEl typically for a given 
direction of k one of the two branches Im{uj)± takes negative values for |k| between 
and a critical value k* at which it passes the zero. The two-dimensional plot are 
given in Figs. El and lAll Other cases are shown in Fig JBll and discussed in Appendix 
[Appendix B| together with the orient at ion- dependence. The critical point k* (a contour 
in 2- dimensions, see Figs. |S] and IA1|) , fixes the minimal unstable wavelength. We 
therefore expect to find unstable solutions associated with the formation and evolution 
of ripples with wavelengths A > A* = 27r//c*. 



5. Ripple wavelength 

Several analytical solutions of EqEl can be found in some special cases which are 
discussed in Section IHl But the study of the surface instabihties can be highly simplified 
if we consider the first order effects when the sputtering flux r] is small. 



5.1. Approximate equation 



In the case of small sputtering fluxes, the branch of Im{uj{k)), with negative values can 
be approximated to: 

Plfc6 + P^k'' + P3(k)fc2 + + P,{k) 



Im{uj) 



V^k^ + 27^/^2 + (vk)2 + 72 



(10) 



with 



When k = 
order obtaining: 



Pi = V[{1 - 0)7/Ci - V{Vi - (1 - 0; 

P2 = {l-(j)h[V{V2-S,)+^IC,]-{ 
P3(k) = -[Pi + (l-0)P2](vk)2 

P4 = 7'[(l-0)si-0I?i] 

P5(k) = (l-0)7(vk)[(Vi-V2)k] 

= |k| is sufficiently small { k <^ /rj 



Ak^ - Be 



1 + 0)7PPi 



(11) 

we can develop EqlTUl at the 4*^^ 

(12) 



with 



A 



^X> + f 2 



r 



v{vi - V2) 



2-fV + V 



2-, 



B = (f)V, + {l- 



si + 



V[Vl - V2j 

7 



(13) 



Here v, vi and V2 are respectivelly the components of v, vi and V2 in the direction 
parallel to k (i.e. they are: v = |v| cos(vk); vi = |vi| cos(vik) and V2 = |v2| cos(v2k)). 
A comparison between this approximate solution and the exact one is given in FigEJ 
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Note that this approximate expression admits a negative minima, hke the one shown 
in FigEl only if the two terms A and B are both larger than zero. Otherwise for A > 
and -B < the negative minimum disappears, whereas for A < the instability moves 
at |k| ^ oo (see Appendix [Appendix B[ for a discussion of the exact case). 



5.2. Solutions 

The expected wavelength of the ripples is associated with the fastest growing mode, 
which corresponds to the value of k at which Im{uj)- reaches its most negative point. 
From Eqll2| the minimum of Im{uj) is at 




Therefore, at the beginning, the roughness will grow exponentially fast as W ^ 
exp{BH / {AA)) with associated ripple-wavelength at: 

A ~ 27r 




5.3. Orientation and Rotation 

The most unstable mode is selected by the position of the deepest minima in the {k^, ky) 
plane. Equation IT^ gives the radial position |k| = A; at which a section of /m(a;(k)) 
along a given direction of k has a minimum. Such a minima corresponds to the value 
/m(ci;(k)) ~ —R^/{4:A) fEqlT^. The azimuthal position a of the absolute minimum is 
along the direction at which this minimum reaches the deepest point (see Figs. El and 
lAlj) . In EqElwe have different sources of orientation- anisotropy. First, the terms Vi 
and 1^2 which depend on the azimuthal angle. Second, the quantities v,Vi,V2 which 
also depend on such an angle. Indeed, they are the components of the respective 
vectors v,fi,f2 along the given direction of k (t> = |v| cos(vk); vi = |vi| cos(vik) and 
V2 = |v2| cos(v2k)). Consequently we expect the formation of ripples along preferential 
directions when by, by and/or dx, dy are anisotropic and/or when the parameters v, 
Vi, V2 are different from zero. Indeed, in Figs. 131 and IXTI we show that the azimuthal 
position of the minimum can rotate by varying the relative weight of the coefficients 
b^ and by. Numerical solutions, reported in Figl?! indicate that ripples form along 
preferential directions and rotations of their orientations can be induced by introducing 
anisotropy in the b^ and by coefficients. Other kinds of ripple-rotations can be induced 
by varying the x-, ^/-components of the vectors v, vi, V2 (Fig JAlj) . An example of 
ripple-rotation induced by the variations of these terms is provided in Fig |A21 

The dependence of the ripple orientation from the incidence angle of the 
ion-sputtering can be easily explained by the associated dependence of the 
adsorption/excavation rates and the drift from the sputtering direction. We also observe 
that in the isotropic case and when the non-scalar quantities v, vi and V2 are set to zero 
the surface instability tends to generate spontaneously pyramidal-like surface-structures 
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Figure 6. Evolution of the surface roughness (log-scale) v.s. time (linear scale) from 
numerical simulations (see Appendix fSppendix D| |. (Arbitrary units) . The insertion is 
the log-log plot of the last part of the evolution (highlighted with the rectangular box in 
the main plot). When the Erlich-Schwoebel barrier is active the dynamical evolution of 
ripples can be described with an exponentially-fast growth at the beginning and then 
a 'saturation' to a slower growth consistent with power-law (linear trend in log-log 
scale). In this figure is also visible a slower behavior at the beginning of the growth- 
process. This transient regime, that often occurs in our numerical simulation, is also 
characterized by an exponential growth but with a longer characteristic time. See 
Appendix [Appendix D| for details on the parameters used in the simulation. 



(see Fig JA2|l . For a further discussion of the orientation-dependence of the instabihty 
modes see Appendix [Appendix B 



5.4- Some special cases 

Let us first observe that, when /Ci, si and (j) are equal to zero, the ripple wavelength, 
given by Eq J15[ coincides with the one found for sand dunes in deserts (see for instance 
[IH]). In our notation the 'reptation length' is /q = v/'j, the 'cut-off length' is 
Ic = {V2 — Vii/v, whereas Vi — V2 is the collective drift velocity of the dunes. The 
approximations usually applied in this context [T71 HH], imply: ^ \/'^ H-, and 
7/c ^ f 1 — f2, giving, from EqlTHl 



. (16) 

Let us now consider the dynamical evolution of a surface under ion sputtering and 
in particular the case when the effect of the Erlich-Schwoebel barrier is not present (as 
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Figure 7. An example of ripples rotation associated with the relative variation of 
the two coefficients bx and by. See Appendix [Appendix D| for details. 



for semiconductors and glasses). In this case, s = 0, Si = and we also expect that 
the drift velocity v and the dispersion constant V are equal to zero or infinitesimally 
small. Indeed, here the current of mobile atoms on the surface is mainly induced by 
the differences in the chemical potential. Under these assumptions, from EqlTHj the 
wavelength of the most unstable ripple is: 



A ~ 27r 



2/C 



(17) 



where we called u = -T>i(f)/{l — (f)), a quantity which plays the role of an effective surface 
tension. Note that Eq. El is the same result as from the Bradley and Harper theory 

j2nil22ll2ZIE3Ell- 

When the Erlich-Schwoebel barriers are active {s,Si ^ 0), effects can be observed 
on the ripple-wavelength at their beginning, which becomes: 



A ~ 27r 



2/C 

u + s 
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6. Exact solutions 



Compact analytical expressions for the values of k at which Im{uj) = (k*) can be 
calculated from Eq|71in some special cases. 

In particular, when = 0, Si = 0, /Ci = and V = 0, we obtain 



where v, vi and V2 are the components of v, vi and V2 in the direction of k*. Whereas 
the value k* is the radial component of k along the closed contour in the {kx, ky) plane 
at which the surface /m(a;(k)) crosses the zero (see Figs. El and Ell- 

On the other hand when, /Ci, Si, Vi and T>2 are equal to zero, we find 



h* - I ~ ^^^^ (on\ 

" - V{v, - V - [l - ct>)v2) ■ ^ ^ 

The effect of the deterministic diffusion induced by the chemical potential can be 
studied from the solution 



^/ :i-0)7/Ci-P(I?i-(l-0M ' ^ ^ 

which holds when v = 0, vi = 0, V2 = 0, si = and P - + (1 - (p)V2 > 0. 
Whereas, when we set to zero vi, V2, T>i and V2, we find 

which implies that the uphill current due to the Erlich-Schwoebel barrier can generate 
instability even when the shape- dependent erosion and recombination terms are inactive. 



7. Exponential growth, stabilization, saturation and critical roughness 

In metals, when the Erlich-Schwoebel barrier is active, there is an important non-linear 
contribution in the current of mobile atoms which becomes sizable when the roughness 
becomes sufficiently large and therefore ((a^V/i)^) ~ 1 (see EqE)) (the average is over the 
surface positions). We observe -numerically- that this term changes the ripple's growth 
dynamics: from exponential to a slower growth consistent with a power-law. Sometime 
a saturation to a constant value is observed. This effect is shown in FigsEl and IHl 
Numerically, all the computed exponents follow in a wide range between (saturation) 
to 1. We observe that this power-law regime is strongly affected by the presence or 
the absence of a noise term in EqI21 In particular the saturation to a fixed profile 
with constant roughness has been observed only in absence of the noise term. In the 
power-law-like regime the ripple-wavelength tends also to grow with time. A theoretical 
justification of this power law regime and the evaluation of the exponent is under current 
investigation. From preliminary studies (described in Appendix [Appendix CD it seems 
apparent that this phase of the growth is characterized by a strong non-linearity with a 
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Figure 8. During the exponential growth the profile shows a typical sinusoidal 
shape (a). Whereas at the saturation/stabilization, the profile tends to eliminate the 
curvature creating more triangular-like shapes (b) . At this stage instabilities might take 
place, with short wave-length profile corrugations which nucleate on the triangular-like 
profile (c). (See text.) 

strong dependence of the ripple shapes (FiglHI) and dynamics (FigslHlandini) from small 
variations of the system parameters or the starting conditions. These preliminary studies 
seem to indicate that this regime is governed by the re-nucleation of new instabilities 
over the rippled surface. Experimentally, power law growth of the roughness and growth 
of characteristic wavelengths were observed in erosive sputtering |3]. On the other hand, 
re-nucleation of instabilities has not been reported yet. 

In semiconductors or glasses, when no Erlich-Schwoebel barrier is present, it is 
physically intuitive that the exponential growth of the surface roughness (which is a 
characteristic of the beginning of the surface instability) cannot continue indefinitely. 
Indeed, from the expression R{r,t) = Rq + Riexp{iut + ikr), which we used to 
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Figure 9. When instabilities re- nucleate on the triangular- like profiles the growth is 
characterized by a power-law regime. (See text.) 

derive EqUHl we can immediately observe that when Ri > Rq = {1 — 0)?7/7, the 
amount of mobile atoms might become negative. Since a negative amount of atoms 
is physically impossible, the process of exponential roughness-growth described above 
must necessarily finish around a critical roughness given by: 

W,^{l-<f))1- . (23) 

7 

This behavior is confirmed by numerical solutions of Eq^ [S3] and it is expected to be 
observable in semiconductors and glasses after sufficiently long times. 

We note that in the present section, these nonlinear effects are only sketched and 
this might lead to some incompleteness and trivialities. The systematic study of the 
nonlinear behavior is a worthwhile exercise which requires further careful investigations 
and it would be the subject of future publications. 

8. Conclusions 

We have shown that the same theoretical approach introduced to describe the formation 
of aeolian sand ripples can be conveniently applied to the study of the formation of 
periodic structures on surfaces under ion sputterning. We have two physically different 
phenomena which both take place at a surface and involve two kind of particles. In 
both cases these particle are excavated, adsorbed and displaced on the surface by 
different physical agents. The equations used to describe the Excavation and Adsorption 
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mechanisms (Eqs. 121121) are the simplest expressions that take into account a dependence 
of these phenomena to the surface shape. Indeed, they are composed by a constant 
term, plus a first derivative, plus a second derivative (i.e. height /orient at ion, slope, and 
curvature). Such terms must appear in any theory that addresses the problem of surface 
excavation and adsorption. On the other hand, the mobility term contains more specific 
terms like the deterministic diffusion or the Erlich-Schwoebel barrier term, which are 
specific of the ion sputtering surface physics. But, also in this case, these terms describe 
a rather simple dependence of the particle mobility form the surface shape. Similar 
terms might be profitably introduced in the description of sand dune dynamics. 

We performed a linear stability analysis for pattern formation under sputtering 
erosion and we obtained general expressions for the ripples wavelength in term of the 
system parameters. It has been shown that in some particular cases such solutions 
coincide with the ones already known in the literature for sand dunes and surface 
instability [201121122111111111 • We have discussed the effect of the Erhch-Schwoebel 
barriers and compared the result with numerical solutions. We pointed out the Erlich- 
Schwoebel barrier can be responsible for a dramatic change in the system dynamics: from 
the exponential growth to a slower growth compatible with a power-law dynamics. The 
occurrence of a critical roughness has been predicted. The effect of ripple-renucleation 
in this regime has been also highlighted. The dependence of the ripple-orientation 
and shapes in function of the system parameters has been discussed and the intriguing 
phenomenon of 'ripple-rotation' has been accounted. 

It should be noted that the main purpose of this paper is to point out a relevant 
example of universality: two processes which have completely different scales present 
a dynamical evolution which obeys to the same geometrical constraints and thus can 
be described by using the same phenomenological model. On the other hand, we must 
observe that the class of solutions of Eq^ is rich and complex - even in the linear 
approximation. Exhaustive, systematic studies of the classes of solutions of this equation 
and their dependence on the set of parameters will be the subject of future studies and 
publications. 
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Appendix A. Fourier transform of the linearized equation 

By substituting Eqs|Hl S IHl and |21 into EqQ and by neglecting the second order terms 
(in Ri and hi), we obtain the following linearized equation: 

^ = 7i?i - [vi - (1 - 0)V2] V/ii - 
[Vi-{l-<f))V2\ V'h 
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dRi 



- 7i?i - vVi?i + VV^Ri + 



dt 



(1-0) (Vi-V2)V/li + 



iVi-V2- si)V^hi- JCiV^h 



(A.l) 



A Fourier analysis of Eq iA.ll leads to 



7^1 - lieu + ik [vi - (1 - 0)v2] - 



fc2 [V^-{l-<p)V2]yh = 
[icj + 7 + ikv + k'^V] Ri-{1- 0) ik(vi - vs) 
p (Di - ©2 - Si) - k^ICi] h = , 



(A.2) 



with Ri and /ii the Fourier components of Ri and hi respectively. This equation is 



a simple linear equation in two variables. It admits a non-trivial solution when the 
determinant of the coefficients is equal to zero. This leads to EqEl 

Appendix B. Orientation dependence of the instability 

There are three distinct mechanisms that can lead to ripple-rotation. 

The first is associated with the anisotropics in the curvature-dependent erosion 
coefficients bx, by and in the curvature-dependent adsorption coefficients dx, dy. These 
anisotropics are induced both by the sputtering orientation angle and by the crystal 
orientation (in non-amorphous surfaces). Example of the instability rotation and ripple 
re-orientation are given in Figs El and [7| 

The second mechanism which can induce ripple-rotation is associated with the 
vectorial terms Vi and V2. In Fig iBll are plotted various solutions of EqEl for the 
branch of Im{uj) which can admit negative values. The plot is v.s. |k| along a given 
azimuthal direction. In Fig lBlb . we show the weakening and disappearance of the 
negative minima in a given direction caused by the increase of the component of vi 
along this direction. In Fig JBlb . we show the re-appearance of the negative minima 
by inceasing the component of V2 along the same direction. Whereas in Fig lBlb . the 
minima is recovered by decreasing simultaneously the components of both vi and V2 
along this direction. In these plots we used: 7 = 1, = 0.5, v = 2, V = 0.1, /Ci = 0.4, 
Si = 7, Vi = 0.7 and V2 = 1. Fig lBlb varies Vi between 4 and 100 at fixed V2 = 1. 
Fig lBlb fixes vi at 100 and increase V2 from 10 to 98. Fig lBlb decreases vi from 100 to 
2 and decreases V2 from 1 to 0.02. 

An analogous example is given in Fig lAll where the position of the minimum rotates 
by approximately 90 degrees by changing the components x, y of the coefficient vi. The 
sensitivity of the instability to the relative weights of the components of the coefficients 
vi, V2 along a given direction is rather evident. 
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(a) (b) 

Figure Al. Surface-plot and contour lines plot for the imaginary part of the 
dispersion relation Im(uj) from the analytical solution of Eq|Sl in the two-dimesional 
case. Figures (a) and (b) show the occurrence of a 90 degrees rotation in the orientation 
of the instability mode associated with the variations of the components in the x, y 
plane of the coefficient vi. To obtain the plots on the left side we set: Vx = 0.1, = 0.1, 
{vi)x = 10, {vi)y = 0.001, V2 0, s = 7, 7 = 1, 7] = 1, 2? = 0.1, , Pi = 0.7, 
1)2 — 1, ICi — 0.4, 4> = 0.2. The right plots use the same parameters with exchanged 
components for vi. 

The third mechanism which influences the orientation of the ripples is associated 
with the drift velocity v. We expect this last mechanism to be more relevant in the case 
of aeolian dunes and less important in ion-sputtered surfaces. 

What we want to stress is that the phase diagram associated with these solutions 
is very rich and non-trivial. Instabilities in a given direction can be trigged on and 
off by changing the relative weights of the quantities b^, by or d^, dy or by varying 
the intensities or the components of the coefficients v,vi,V2. These changes could be 
associated with variations of the angular orientation of the ion sputtering. It is beyond 
the propose of the present paper to discuss in detail these aspects, however we want to 
make clear that the present theory can account these phenomena. 
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Figure A2. An example of ripples rotation (left and right figures) associated with the 
variation of the components of v, vi, V2 in the {x, y) plane. Different profile instability 
(from modulated rippled-structure to pyramidal-like surface-structure) might take 
place when the parameters v, vi and V2 go to zero (central figure). 




Figure Bl. The imaginary part of the dispersion relation Im{Lo) v.s. |k| along a 
given azimuthal direction for various values of the components of v, vi , V2 . Instabilities 
become deeper or rise and even disappear depending on the weights of the components 
of v,vi,V2 (see text). 
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Appendix C. Saturation/Stabilization 

We observed numerically that the behavior of the surface instability in the 
saturation/stabilization regime described in Section [7| is very sensitive to the presence 
of a noise term added to Eq|2l (see main text). In order to investigate this effect we 
performed several simulations by varying the amplitude (a) of an additive Gaussian noise 
with average zero and variance . Such a parameter a was varied from zero to an upper 
limit which has been chosen several order of magnitude smaller of the critical roughening 
Wc- We observe that the stabilization/saturation dynamics is affected by the amplitude 
of the additive noise. In particular, in some region of the parameters we observed 
that the saturation to a constant roughening can be achieved only without any noise 
contribution, whereas the exponents in the power-law growth appear to be dependent 
on the amplitude of the noise. But it must be noted that the observed behavior is 
extremely complex and shows a strong dependence on the starting configuration, on the 
history (time steps, random generated noise etc) and on the system parameters. After 
all, this is not surprising in a non-liner system. 

Let us here briefly describe the kind of (non-linear) dynamics that we observe 
numerically in the one-dimensional case. After the first exponential growth when 
the Erlich-Schweoebel barrier is active, the system reaches a saturation/stabilization 
regime where the roughening grows slower than exponential or stays constant. In the 
exponential-growth stage we observe that the surface-profile has a typical sinusoidal- 
like shape (see FiglHK)- After this stage, at the saturation/stabilization regime the 
surface tends to assume a triangular like profile flattening in this way its curvature and 
concentrating it at the vertices (Fig|H|D). In this regime the growth of the surface-ripples 
might end or -vice versa- it might re-start from small instabilities (infinitesimal noise) 
generated over the triangular shapes (FiglSl:;). In this second case, the re-nucleated 
ripples follow a similar history of the ones on which they nucleate: they start with and 
exponential growth and then reach a saturation until another re-nucleation occurs. In 
this phase, interactions among ripples, nucleated in different part of the surface, and 
interaction between the newly nucleated ripples and the pre-existent surface instabilities 
plays a major role. Typically, the overall dynamics appears to be consistent with a power 
law growth associated with a change in the ripple wavelength, but a finer analysis of 
this process shows a complex non- uniform growth (see FigEl). 

Figures|HK,b were generated by setting: v = 0.1, V = 0.2, K. = 3, s = 0.6, = 10~^, 
7] = 0.05, 7 = 0.03, a = 1, b = 5, c = 0.1, d = 0.5, a = 10^, with no additive noise. 
Whereas Figs. |Hfc,inihave the same set of parameters but with an infinitesimal additive 
Gaussian noise C, (EqEj) with zero mean and standard deviation equal to 10~^^. 

Appendix D. Numerical Solutions 

The numerical solutions of EqQ presented in this paper and in particular the ones shown 
in FigsIS El and IHl have been performed as follows. We considered a one-dimensional 
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flat substrate {h{x, 0) = hO) of lengtli L, witli periodic boundary conditions. An 
infinitesimal quantity of mobile atoms were added randomly to the substrate (with 
< R{x, 0) < L/N10~^^). We then computed the profile-evolution using EqlUwith the 
derivative substituted with finite differences. To this purpose, the substrate has been 
divided into discrete points. The -adimensional- time indicated in FiglHl and El is the 
number of numerical steps. The height is in unit of L/N and the roughness is defined 
as w{t,L) = {[h{x',t) - {h{x,t))J'^)l('^ (see, for instance, 

Several computations with a number of points equal to = 1000, 2000 and 3000 
(the one presented here have A^ = 3000) have been performed to verify the effect 
of boundary and discretization. Moreover, simulations with no periodic boundary 
conditions and with the sputtering term (EqI21) applied only to a central mask, have also 
been performed obtaining very similar results. The robustness of the present approach 
has been verified varying the parameters, the time steps, the initial roughness of the 
substrate, etc. Comparable results have been always found but, we must stress that, 
under some conditions, numerical instabilities (in particular surface-deformations with 
A ~ L/N) can be trigged on depending on the protocol utilized. 

The simulation result shown in FigsEl IHl and El use: v = 0.1, V = 0.2, /C = 3, 
s = 0.6, = 10-^ ri = 0.05, 7 = 0.03, a = 1, b = 5, c = 0.1, d = 0.5, a = 10^ The 
simulation time was 30000, steps. 

The numerical solutions for the two-dimensional case have been calculated by 

following the same protocol as in the one-dimensional case. In two dimensions the 

substrate is a square L x L surface with periodic boundary conditions which has been 

subdivided into a discrete grid of 150 x 150 points. The simulation result shown in FiglSl 

uses: V = 0, 1) = 0, /C = 0.5, s = 0, (f) = 0.2, r] = 0.05, 7 = 0.03, a = 0, 6^ = 2, 6^ = 2 

c = 0, dx = 0.3, dy = 0.3, a = 0. = 2, by = 5. Whereas the simulation in Fig|71uses 

the same parameters except for bx = 2, by = 5 (up) and b^ = 5, by = 2 (down). The 

simulation time was 5000, steps. 
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